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As a neutron star (NS) is tidally disrupted by a black hole (BH) companion at the end of a BH- 
NS binary inspiral, its magnetic fields will be stretched and amplified. If sufficiently strong, these 
magnetic fields may impact the gravitational waveforms, merger evolution and mass of the remnant 
disk. Formation of highly-collimated magnetic field lines in the disk+spinning BH remnant may 
launch relativistic jets, providing the engine for a short-hard GRB. We analyze this scenario through 
fully general relativistic, magnetohydrodynamic (GRMHD) BHNS simulations from inspiral through 
merger and disk formation. Different initial magnetic field configurations and strengths are chosen 
for the NS interior for both nonspinning and moderately spinning (aBH/MBH=0.75) BHs aligned 
with the orbital angular momentum. Only strong interior (Bmax 10^^ G) initial magnetic fields in 
the NS significantly infiuence merger dynamics, enhancing the remnant disk mass by 100% and 40% 
in the nonspinning and spinning BH cases, respectively. However, detecting the imprint of even a 
strong magnetic field may be challenging for Advanced LIGO. Though there is no evidence of mass 
outfiows or magnetic field collimation during the preliminary simulations we have performed, higher 
resolution, coupled with longer disk evolutions and different initial magnetic field configurations, may 
be required to definitively assess the possibility of BHNS binaries as short-hard GRB progenitors. 

PACS numbers: 04.25.D-,04.25.dk,04.30.-w 



I. INTRODUCTION 

With the first direct detection of gravitational waves 
(GWs) expected in the next few years, numerical rela- 
tivity simulations will be crucial for distinguishing dif- 
ferent GW sources from one another. Mergers of black 
hole-neutron star (BHNS) binaries are among the most 
promising sources of gravitational waves detectable hy 
ground-based laser interferometers like LIGO [H, 
VIRGO [31,14], GEO LCGT [6], and AIGO [7^, as well 
as by the proposed space-based LISA-like interferome- 
ters and DECIGO [9] and third-generation ground- 
based detectors such as the Einstein telescope [1Q|, [Tlj . 
Analysis of gravitational waveforms from BHNS mergers 
may spark new insights into the behavior of matter at 
nuclear densities. 

Theoretical models indicate that a neutron star- 
neutron star (NSNS) [l^-^ or BHNS [is", '20'-'26] binary 
merger may result in a hot, massive disk around a black 
hole, whose temperatures and densities could be suffi- 
cient to trigger a short-hard gamma-ray burst (SGRB). 
Indeed, SGRBs have been associated with galaxies with 
extremely low star formation rates (see j27j and refer- 
ences therein for a review), indicating that the source is 
likely to involve an evolved population, rather than main 
sequence stars. The number of detectable BHNS merg- 
ers in the observable universe is still an open question, 
due to uncertainties in population synthesis calculations. 
The estimated event rate of BHNS mergers observable by 
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an Advanced LIGO detector typically falls in the range 
7^ - 0.2-300 yr-i [28, 29]. 

Motivated by the significance of BHNS binaries both as 
detectable GW sources and SGRB candidates, many sim- 
ulations of BHNS systems have been performed in past 
years in a Newtonian or post-Newtonian framework (see, 
e.g., jSOMsj) and in conformally-fiat relativistic gravita- 
tion |20|, [21| . Recently, several groups have performed dy- 
namical simulations of BHNS binary inspirals and merg- 
ers in full GR [liM llil^- 

Over the past few years, we have studied BHNS merg- 
ers beginning with the construction of quasiequilibrium 
circular orbit initial data ^48n52i| a nd following up with 
full GR dynamical simulations [lH, |26| . Our GR simula- 
tions, and those by other groups, suggest that for initially 
nonspinning BHs, the remnant disk mass is substantial 
for q = Mbh/^ns ^ 3 and tends to increase with de- 
creasing q. For a fixed q^ the disk mass also increases 
for smaller NS compaction and for more rapidly spin- 
ning BHs aligned with the orbital angular momentum of 
the binary. For sufficiently high spins, small mass ra- 
tios, and/or lower NS compactions, a substantial disk 
can form following the merger, favoring BHNS mergers 
as plausible central engines for SGRBs. However, these 
simulations have yet to account for magnetic field effects 
- a crucial component in many SGRB models involving 
a disk around a spinning BH (see, e.g. [53-55]). 

For NS surface field strengths B < lO^^G, magnetic 
fields are unlikely to affect the dynamics of the BHNS 
inspiral and merger [56]. This was shown to be the case 
for NSNS binary inspirals in (57|-[6Q|. Despite this con- 
clusion, magnetic fields may significantly influence the 
post-merger dynamics, as the fields are likely to be am- 
plified during and after merger. Magnetic fields could 
stir turbulence in the remnant disk, resulting in angular 
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momentum transport and accretion onto the BH. They 
could also lead to matter outflow and jets along the BH 
remnant spin axis 1611. another ingredient required by 
most SGRB models |53M55|. 

In this paper, we present a new set of fully relativis- 
tic BHNS simulations that probes how magnetic fields 
infiuence the dynamics and outcome of the merger us- 
ing our new adaptive mesh refinement (AMR) GRMHD 
code [H, [63|. Fixing the BH:NS mass ratio at = 3, 
we consider the cases where the BH possesses no spin 
(A cases) and moderate spin aBn/^BH = 0.75 aligned 
with the orbital angular momentum (B cases). Since the 
internal magnetic field strength and configuration in a 
NS is not known, we vary the strength and geometry of 
the internal fields to study their effects. We find that 
for low and moderate field strengths ^ lO^^G, magnetic 
fields do not significantly alter the inspiral and merger 
dynamics, which is consistent with the result reported 
in [42). Here, by low and moderate magnetic fields we 
mean those with field strengths small or moderate when 
compared to the virial value of about lO^^G; well below 
this value the field is dynamically unimportant. However, 
when the central field strength approaches ~ lO^^G, cor- 
responding to a magnetic to gas pressure ratio of ~ 0.5%, 
the merger dynamics and remnant disk mass are affected 
significantly. Yet even with such strong internal mag- 
netic fields, the emitted GWs are not appreciably differ- 
ent from the unmagnetized case, at least for the prelim- 
inary set of models considered here. In the late inspiral 
and merger phases, tidal deformation and disruption of 
the NS play key roles in distinguishing GWs from BHNSs 
and BHBHs. 

During merger, most of the magnetized NS matter is 
captured by the BH. Only when the NS interior is seeded 
with strong magnetic fields (5max ^ lO^^G, near the 
center of the NS) is a significant impact on the dynamics 
observed, resulting in a disk that has up to twice the 
rest mass as the corresponding unmagnetized case. In 
all cases the disk accretion rate onto the BH decreases 
with time immediately after merger, before settling down 
to a quasistationary state. Most of the magnetic field 
lines are tightly wound within the remnant disk, and no 
evidence of magnetic field collimation around the final 
spinning BH is observed by the time we terminate our 
simulations. The remnant disk is hot (T ^ IMeV) and 
massive (Mdisk ^ O.O2M0 and ~ O.IM© for cases A and 
B, respectively). 

The magnetic fields threading the remnant disk may 
be amplified and tangled on a longer timescale than we 
simulate, stirring up MHD turbulence. Based on extrap- 
olation of the accretion rates near the end of our simula- 
tions, the lifetime of the disk is roughly O.3(Mo/1.4M0)s. 
While there is no evidence of outfiows during these pre- 
liminary simulations, longer disk evolutions, higher res- 
olution and different B-field geometries may be required 
to definitively assess the possibility of BHNS binaries as 
short-hard GRB progenitors. 

The following sections are organized as follows. Sec- 



tions ini and Uni review the basic equations, including 
our initial data, gauge conditions, matter evolution equa- 
tions, and diagnostics, as well as their implementation in 
our GRMHD code. Section [IVl presents the results of our 
magnetized BHNS merger simulations. Finally, we sum- 
marize our findings and comment on future directions in 
Sec. El 



II. BASIC EQUATIONS 

This section introduces our notation, summarizes our 
method, and points out the latest changes to our numer- 
ical technique as summarized in [25|, [26|, [62|, |63|. Ge- 
ometrized units {G = c = 1) are adopted, except when 
stated otherwise. Greek indices denote all four space- 
time dimensions (0, 1, 2, and 3), and Latin indices label 
spatial parts only (1, 2, and 3). 

We use the 3+1 formulation of general relativity and 
decompose the metric into the following form: 

ds^ = -a^dt^ + -fij{dx' + p'dt){dx^ + p^dt) . (1) 

The fundamental variables for metric evolution are the 
spatial three- metric 7^^ and extrinsic curvature Kij. 
We adopt the Baumgarte-Shapiro-Shibata-Nakamura 
(BSSN) formalism [64, 65] in which the evolution vari- 
ables are the conformal exponent (j) = ln(7)/12, the con- 
formal 3-metric 7^^ = e~^^jij^ three auxiliary functions 
= — 7*-^ ,j, the trace of the extrinsic curvature and 
the trace-free part of the conformal extrinsic curvature 
Aij = e-^'^{Kij - -fijK/'i). Here, 7 = det(7i^). The fuh 
spacetime metric g^jj is related to the three-metric 7^^^ 
by 7/izy = Qjjiv + where the future-directed, timelike 

unit vector normal to the time slice can be written in 
terms of the lapse a and shift as = a~^{l^—j3^). 
Evolution equations for these BSSN variables are given 
by Eqs. (9)-(13) in We adopt standard puncture 

gauge conditions: an advective "1+log" slicing condi- 
tion for the lapse and a "F-freezing" condition for the 
shift [66] . The evolution equations for a and are given 
by Eqs. (2)-(4) in [26], with the r] parameter set to 2.2/M 
for the initially nonspinning BH cases and 3.3/M for the 
spinning BH cases, where M is the ADM mass of the 
BHNS binary. We add a fifth-order Kreiss-Oliger dissi- 
pation term to all evolved BSSN, lapse and shift variables 
to reduce high-frequency numerical noise associated with 
AMR refinement interfaces (see [67] for a review and ref- 
erences). 

The fundamental MHD variables are the rest-mass 
density po, specific internal energy e, pressure P, four- 
velocity u^^ and magnetic field = nyF""^^. Here F*/^^ 
is the dual of the Faraday tensor F^^ . Note that is 
purely spatial {B^ = —n^B^ /a = 0). We adopt a F-law 
equation of state (EOS) P = (F - l)poe with F = 2, 
which reduces to an n = 1 polytropic law for the initial 
(cold) neutron star matter. The stress-energy tensor is 
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given by 

T^u = {poh + h^)u^Uy + + g^y - h^hy , (2) 
where = 1 + e + P/po is the specific enthalpy and 

is the magnetic field measured in fiuid's comoving frame, 
modulo a factor of l/>/47r. Here P^^ = g^^ + u^u^ and 
^2 _ 5^5^ xhe comoving magnetic energy density is 
u^u^T>^\, = 62/2, where T^l, = b'^u'^u^ + (672)5'"^ - 
6^6^ is the stress-energy tensor associated with the mag- 
netic field. In the ideal MHD limit, in which the plasma is 
assumed to have perfect conductivity, the Faraday tensor 
can be written as F^^ = u^e^^^^bs. 

In the standard numerical implementation of the MHD 
equations using a conservative scheme, it is useful to in- 
troduce the "conservative" variables p*, 5^, f and P*. 
They are defined as 

P* -V^Pon^u^ , (4) 
^ -VlT.yn^r, , (5) 
r = ^T^un^n^-p, , (6) 

^ V7P^ (T) 

The evolution equations for p*. Si and f can be derived 
from the conservation of rest mass V ^{p^u^) — and the 
conservation of energy- momentum V^T^^ = 0, giving 
rise to Eqs. (27)-(30) in [62|. 

In the ideal MHD limit, the Maxwell equation 
VzyP*^^ _ Q yiei(jg the magnetic constraint d^B^ — 
and induction equation 9^P* + dj{v^ B'^ — v'^&) = 0. As 
shown in [68| and [62], these equations can be rewrit- 
ten by introducing the electromagnetic 4-vector poten- 
tial Afj^ = + A^, with n^A^ = 0. The magnetic 
constraint and induction equations become 

P^ = e'^^djAk , (8) 
dtA^ = e^jkv'B^ - di{a^ - (3^ Aj) , (9) 

where e*-^^ = n^e^*-^^ is the 3-dimensional Levi-Civita 
tensor. In we evolve the vector potential and choose 
the algebraic EM gauge ^ = P^Aj/a = —n^Aj. We 
have found that for BHNS simulations, we can achieve 
better numerical results by imposing the Lorenz gauge 
^fiA^ = 0, which gives the evolution equation 

dtiV^^) + dj{a^A^ - = (10) 

for <l> in place of the algebraic gauge condition [63] . Our 
numerical implementation of Eqs. (|8j) and (|9j) guarantees 
numerically identical P* (see [62]) regardless of EM gauge 
in simulations with a uniform-resolution grid. However, 
interpolations performed on Ai at refinement boundaries 
on AMR grids will modify A^, resulting in different P* 



near these boundaries. We have shown that in the alge- 
braic gauge ^ = —n^Aj, there exists a zero-speed mode, 
which in BHNS simulations manifests itself as a trail of 
nonzero Ai left behind the orbiting NS [63j. When AMR 
refinement boxes tracking the motion of the NS cross this 
"trail", spurious, strong magnetic fields appear on the 
refinement boundaries. On the other hand, the Lorenz 
gauge exhibits no zero-speed modes. As a result, the be- 
havior of the P* fields on refinement boundaries is drasti- 
cally improved [63] . We therefore adopt the Lorenz gauge 
for all simulations in this paper. 



III. NUMERICAL METHODS 
A. Initial data 

Our initial data are constructed by solving Einstein's 
constraint equations in the conformal thin-sandwich 
(GTS) formalism, which allows us to impose an approx- 
imate helical Killing vector by setting the time deriva- 
tives of the conformally related metric jij to zero in the 
frame corotating with the binary. We model the NS as 
an irrotational n = 1 polytrope, and impose the black 
hole equilibrium boundary conditions of Cook and Pfeif- 
fer [69] on the black hole horizon. The CTS initial data 
correspond to a binary in circular quasiequilibrium with 
a separation chosen to be outside the tidal disruption ra- 
dius. Details of this method can be found in [26|,|52|. The 
initial data used in this paper are the same as case A (for 
an initially nonspinning BH) and case B (for an initial BH 
spin Jbh/M|h = described in [26]. 

The initial data are calculated using the Lorene spec- 
tral methods numerical libraries [70] . The excised BH re- 
gion is filled with constraint- violating initial data, using 
the "smooth junk" technique we developed and validated 
in [71] (see also 0,[73|)- In particular, we extrapolate all 
initial data quantities from the BH exterior into the inte- 
rior with a 7*^ order polynomial, using a uniform stencil 
spacing of Ar ~ 0.3rAH- 

All of the NSs considered in this paper have a com- 
paction of C = Mns/Pns = 0.145, where Mns is the 
ADM mass and Pns is the (circumferential) radius of the 
NS in isolation. Since we model the NS with an n = 1 
(F = 2) polytropic EOS, the rest mass of the NS, Mq, 
scales with the polytropic constant k, as Mq cx k,^^'^. For 
a NS with compaction C = 0.145, we find the ADM mass 
for the isolated NS to be Mns = 1.3OM0(Mo/1.4M0), 
with an isotropic radius Piso = 11.2km(Mo/1.4M0) 
and circumferential (Schwarzschild) radius of Pns = 
13.2km(Mo/1.4M0). The maximum rest-mass density 
of this NS is po,max = 9 X lO^^^g cm-3(1.4M0/Mo)^ 

We add a small, poloidal magnetic field via the vector 
potential of the form 

A, = (-y^5\ + '^^5y^ (11) 
A^ = Abwl max(P - Peut, 0)"^ (12) 
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FIG. 1. Magnetic field lines (white) at the time when they 
are seeded into the NS for case A4. The NS density profile 
is shown such that darker colors indicate higher densities. 
The BH apparent horizon is shown in black, and the BH-NS 
coordinate separation is 6.12M 

where (xc^^O) is the coordinate location of the center 
of mass of the NS, = {x — Xc)^ + (^ — Vc)^ ^ and A5, 
Tirp and Pcut are free parameters. The cutoff pressure pa- 
rameter Pcut confines the B-field inside the neutron star 
to reside within P > Pcut- The parameter 715 determines 
the degree of central condensation of the magnetic field. 
Similar profiles of initial magnetic fields are also used 
in numerical simulations of magnetized accretion disks 
(see, e.g. [7l,l75j ) and magnetized compact binaries (see, 
e.g. |42|, l76l - [78| ). We set Pcut to be 4% of the maximum 
pressure and 77,5=1 or 2. The parameter A\) controls the 
strength of the initial magnetic field, which can be char- 
acterized by the maximum magnetic field inside the NS, 
as well as the magnetic energy hA defined as 

M = ^ n^n.T^^dV, (13) 

where dV = ip^d^x is the proper volume element on a 
t = constant spatial slice. A4 is the EM energy measured 
by a normal observer. 

Since the magnetic field is expected to remain frozen 
into the NS during the inspiral phase, we add the fields 
immediately before tidal disruption to minimize numer- 
ical error. Magnetic fields added to the NS at t = 
maintain their original profile within the star for much 
of the first orbit {t ^ milliseconds). Table |I] summarizes 
the initial data used in our simulations. Figure [1] shows 
the magnetic field configuration for one of the cases (A4) 
at t = 448. 5M, the time at which the NS is seeded with 
magnetic fields. The seed magnetic fields are too weak 
to significantly perturb the quasiequilibrium NS, leading 
to virtually no change in gravitational field constraint 
violations. 



B. Evolution of the metric and MHD 

We evolve the BSSN equations with fourth-order ac- 
curate, centered finite-differencing stencils, except on 
shift advection terms, where we use fourth-order accurate 



upwind stencils. We apply Sommerfeld outgoing wave 
boundary conditions to all BSSN fields. Our code is em- 
bedded in the Cactus parallelization framework [80[, and 
our fourth-order Runge-Kutta timestepping is managed 
by the MoL (Method of Lines) thorn, with a Courant- 
Friedrichs-Lewy (CFL) factor set between 0.0625 (in the 
coarsest refinement level) and 0.5 (in the innermost 4 re- 
finement levels) in all simulations. We decrease the CFL 
factor in the coarse refinement levels so that we can use 
a larger value for the parameter 77 in the shift equation 
[Eq. (4) in [26]]. We use the Carpet ^ infrastructure to 
implement the moving-box adaptive mesh refinement. In 
all AMR simulations presented here, we use second-order 
temporal prolongation, coupled with fifth-order spatial 
prolongation. The apparent horizon (AH) of the BH is 
computed with the AHFinderDirect Cactus thorn [82|- 
The GRMHD equations are evolved by a high- 
resolution shock-capturing (HRSC) technique j83| that 
employs PPM [84] coupled to the Harten, Lax, and 
van Leer (HLL) approximate Riemann solver [85|. The 
adopted MHD scheme is second-order accurate for 
smooth flows, and first-order accurate when discontinu- 
ities (e.g. shocks) arise. To stabilize our scheme in regions 
where there is no matter, we maintain a tenuous atmo- 
sphere on our grid, with a density floor patm set equal 
to 10~^^ times the initial maximum density on our grid. 
The initial atmospheric pressure Patm is set equal to the 
cold polytropic value Patm = ^Patm- Throughout the 
evolution, we impose limits on the atmospheric pressure 
to prevent spurious heating and negative values of the 
internal energy e due to numerical errors. Specifically, 
we require Pmin 

< P < 

Pmax5 where Pmax — lOtvpQ and 
Pmin = Whenever P exceeds Pmax or drops below 

Pmin, we reset P to Pmax or Pmin, respectively. Applying 
these limits everywhere on our grid would artificially sap 
the angular momentum in the tidally disrupted NS, al- 
lowing matter to fall spuriously into the BH and thereby 
suppressing disk formation [26]. To effectively eliminate 
this spurious angular momentum loss, we impose these 
pressure limits only in regions where the rest-mass den- 
sity remains very low {po < lOOpatm) or deep inside the 
AH, where i/j^ > ip^^^ as in [26]. Here i/j = and we set 
ip^Yir between 10 and 30. 



C. Evolution of magnetic field 

We evolve the magnetic induction equation via the 4- 
vector potential using Eqs. ([9]) and (p!Q|) . We stagger 
the Ai and P* as in [62]. We store ^ on a staggered 
grid (i+,j+,/c+) [all the other hydrodynamic, BSSN, 
lapse and shift variables are stored at (z,j, /c)], where 
z+ = i + 1/2 and similarly for j+ and . We treat the 
term —dj{/3^-^^) in Eq. (p!Q|) using a second-order up- 
wind scheme. We evolve Eq. (j9j) using the finite- volume 
equations similar to Eqs. (63)-(65) in [62[, modified to 
take into account the second term in Eq. ([9]), which does 
not vanish in the Lorenz gauge. The detailed implemen- 
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TABLE I. Initial data for the BHNS simulations. Here asH /Mbh is the BH spin, Q is the orbital angular frequency, Bmax is 
the maximum value of magnetic field inside the NS, assuming the rest mass of the NS is IAMq] A4 is the energy of magnetic 
field defined in Eq. ([13]) at the time when B-field is added; W is the gravitational potential energy of the NS in isolation defined 
in Eq. (65) of 79]. All BHNS systems considered here have the BH:NS mass ratio of 3:1. 
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0.0 


0.0330 


















Al 
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0.04 


2 
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1.3 X 10^^ 
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X 10" 


-5 


0.04 


1 




A3 






1.4 X 10^^ 


1.2 


X 10" 


-5 


0.001 


1 




A4 






9.7 X 10^^ 


5.9 


X 10" 


-4 


0.001 


1 




BO 
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0.0328 


















Bla 






1.3 X 10^^ 


3.1 


X 10" 


-6 


0.04 
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633.6M 
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1.2 X 10^^ 


2.9 


X 10" 


-6 


0.04 
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752.8M 


B2 






1.3 X 10^^ 


1.1 


X 10" 


-5 


0.04 


1 




B3 






1.4 X 10^^ 


1.2 


X 10" 


-5 


0.001 


1 




B4 






9.7 X 10^^ 


5.9 


X 10" 


-4 


0.001 


1 





tation is described in (63| . 

The particular staggering of the Ai and $ variables 
coupled with the particular implementation of the HRSC 
scheme are designed to ensure that the resulting B field 
obtained by taking the curl operator on [Eqs. (60)- 
(62) in [62j] is numerically identical to the standard con- 
strained transport scheme based on a staggered algo- 
rithm [865. We have carefully designed an algorithm for 
the extra term in Eq. (|9j) so that the additional terms 
in the A^ evolution equations cancel exactly after taking 
the curl operator. The resulting numerical values of B^ 
are thus gauge- invariant in unigrid simulations. We have 
confirmed numerically that this is indeed the case. How- 
ever, in simulations with an AMR grid, since we perform 
interpolations on Ai between refinement levels, values of 
Ai are not the same in different EM gauges at the refine- 
ment boundaries. The resulting B field at the refinement 
boundaries is also different in general but should converge 
to a unique, true solution with increasing resolution in 
any gauge. 

As in other numerical relativity simulations, some 
gauges are better behaved than others. We have demon- 
strated in [63] that the Lorenz gauge is superior to the 
algebraic gauge in magnetized BHNS simulations. We 
therefore adopt the Lorenz gauge in all of the magne- 
tized BHNS simulations presented here. 



D. Recovery of primitive variables 

At each timestep, we need to recover the "primitive 
variables" po, and from the "conservative" variables 
p*, f, and Si. We perform the inversion by numerically 
solving two nonlinear equations via the Newton- Raphson 
method as described in [87], using the code developed by 
Noble et al [88]. 

Sometimes the "conservative" variables may assume 



values which are out of physical range, resulting in un- 
physical primitive variables after inversion (e.g. negative 
pressure or even complex solutions). This usually hap- 
pens in the low-density "atmosphere" or deep inside the 
BH interior where high-accuracy evolution is difficult to 
maintain. Various techniques have been suggested to 
handle the inversion failure (see, e.g. [89]). Our approach 
is mainly to impose constraints on the conservative vari- 
ables to reduce the inversion failure. 

One reason for the inversion failure comes from 7^^ 
losing positive-definiteness during tfie BSSN evolution 
due to numerical inaccuracy, wfiicfi occurs only near tfie 
"puncture" deep inside tfie BH. Before performing ttie 
inversion, we cfieck if 7^^ is positive-definite by finding 
its eigenvalues. If 7^^ is not positive definite, we reset 
7ij i^^fij during tfie inversion, where fij is the 3D flat 
metric tensor. 

In the absence of magnetic fields, the inversion failure 
can be avoided completely by enforcing the constraints 
(see [25] and Appendix [A|) 

= -f'^SiSj < f{f + 2/)*), and (14) 
f > , (15) 

which are the necessary and sufficient conditions for the 
inversion to produce the primitive variables in the phys- 
ical range for the T-law EOS with 1 < T < 2 (see Ap- 
pendix [Aj. We enforce these constraints in regions where 
there are no magnetic fields. When the second condition 
is not met, we reset f = fatm = 10"^^fomax, where fomax 
is the maximum value of f initially. When the first con- 
dition is violated we rescale Si so that its new magnitude 
is =f(f + 2p*). 

In the presence of magnetic fields, no simple analogous 
formulae are available. However, one can prove that (see 
Appendix [A]) 

f > V"^^ (16) 
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for any value of primitive variables in the physical range. 
We therefore impose inequality (p!6|) everywhere: if it 
is violated, we reset f = fatm + V^~^^^/87r. However, 
this does not guarantee that the inversion always pro- 
duce physically acceptable primitive variables. If failures 
occur outside the BH after imposing (fT6]) . we apply a fix, 
which consists of replacing the energy equation ([6j) by the 
cold EOS, P = Pcoid(po) = i^Po when solving the system 
of equations, where is the polytropic constant. One 
can show that this procedure always produces physically 
valid primitive variables (see Sec. IA4p . However, we find 
that this fix gives rise to discontinuous data deep inside 
the BH and these data will eventually propagate out of 
the BH horizon. Note that the success of the "smooth 
junk" technique fH] requires that the constraint violat- 
ing initial data filling the interior of the BH horizon be 
smooth. If discontinuous data are used, then information 
can leak out of the BH horizon. To avoid information 
leakage outside the BH, we apply the following condi- 
tions deep inside the BH: 



and 



fm>0, 



(17) 



(18) 



where B' = B'/V^, B'^ = 
the quartic equation 



jijB^B^, and Wm satisfies 



-{B^S^)\B^ 



2Wm)=0. 



It can be shown that the inequalities (fTT]) and (fT8|) are 
sufficient (but not necessary) conditions for the inver- 
sion to yield physically valid primitive variables (see Ap- 
pendix [A]) • We therefore only use them deep inside the 
BH where t/j^ > V^f^r- choose the parameter t/j^^^ be- 
tween 10 and 30. Since the inequalities (pT|) and (fT5|) are 
sufficient but not necessary conditions, we do not impose 
them strictly, but adopt the procedures described at the 
end of Sec. lA 21 We find that this technique gives rise to 
smoother data in the BH interior preventing contraint- 
violating information from leaking out of the BH horizon. 



E. Diagnostics 

During the evolution, we monitor the Hamiltonian 
and momentum constraints calculated by Eqs. (40)-(43) 
of [2^. We also monitor the interior mass Mint and {z- 
component of) the interior angular momentum Jint of the 
system contained in the simulation domain. These quan- 
tities are defined as integrals over the surface of the outer 
boundary dV) of the computational domain: 

Jint = ^e1,- / a;^ (ifr " 5TK)d^m, (20) 
on Jdv 



(19) 



where is the fiat-space Levi-Civita tensor. As pointed 
out in [26| , the integrals can be evaluated more accurately 
by alternative expressions via Gauss's law (67j : 



Mi, 



X 



1 



V \ l^lT ' 



167r ^^'^ 



(21) 



7 — ——^ ^ 



(22) 



where 5* is a surface surrounding the BH horizon, V 
is the volume between S and the outer boundary, p = 
n^rijyT^^, and R is the Ricci scalar associated with the 
conformal 3-metric jij. If our outer boundary were ex- 
tended to spatial infinity, these integrals would yield the 
ADM mass and angular momentum of the system and 
would be constant in time. While our outer boundary 
dV resides in the nearly Minkowski asymptotic regime, 
it is at a finite distance from the BHNS system. Thus 
the integrals are only approximately equal to the ADM 
M and J at t = and decreases with time, due to GWs 
carrying away energy and angular momentum through 
dV. 

When hydrodynamic matter is evolved on a fixed uni- 
form grid, our hydrodynamic scheme guarantees that the 
rest mass Mq is conserved to machine roundoff error. 
This strict conservation is no longer maintained in an 
AMR grid, where spatial and temporal prolongation is 
performed at the refinement boundaries. Hence we also 
monitor the rest mass 



Mo 



/ 



(23) 



during the evolution. Rest-mass conservation is also vio- 
lated whenever po is reset to the atmosphere value. This 
usually happens only in the very low-density atmosphere 
or deep inside the AH where high accuracy is difficult to 
maintain. 

We measure the thermal energy generated by shocks 
via the entropy parameter K = P/Pcoid, where Pcoid = 
hzpQ is the pressure associated with the cold EOS. The 
specific internal energy can be decomposed into a cold 
part and a thermal part: e = ecoid + ^th with 



Ccold 



cold 



d{l/po) 



r-i 



r-1 

Hence the relationship between K and eth is 

1 P Hi 



(24) 



eth = e — ecoid 



r-lpo 
= {K - l)ecoid • 



1 



(25) 
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For shock-heated gas, we always have K > 1 (see Ap- 
pendix B of [26]. 

Finahy, we monitor the mass and spin of the BH during 
the evolution. They are computed using the isolated and 
dynamical horizon formalism [90], with the approximate 
axial Killing vector on the horizon computed as in [91]. 

F. Gravitational wave extraction 

Gravitational waves are extracted using the Newman- 
Penrose Weyl scalar ?/'4 at various extraction radii be- 
tween 50M and 130M. We decompose 7/^4 into s = —2 
spin-weighted spherical harmonics up to and including 
/ = 4 modes. At each extraction radius, the retarded time 
is computed using the technique described in Sec. IIB 
of [92I to reduce the near-field effect. The wavetrain /i+ 
and hx for each mode are computed by integrating the 
corresponding mode of tp^ twice with time usingthe fixed 
frequency integration technique described in |93| . 

We compute the radiated energy AEqw^ z-component 
of angular momentum A Jew and linear momentum 
APq^^ using expressions equivalent to Eqs. (33), (39), 
(40) and (49) of [9^. To check the violation of energy and 
angular momentum conservation, we monitor the quan- 
tities 

5E=[M- M,^,{t) - AEGw{t)]/M , (26) 
5J=[J- Jint(t) - AJGw{t)]/J , (27) 

where J is the ADM angular momentum of the initial 
binary, Mint(^) and Jint{t) are the interior mass and an- 
gular momentum of the system at time t as calculated by 
Eqs. (j2l]) and (|22]). 

IV. RESULTS 

We have performed magnetized simulations of BHNS 
binaries with BH:NS mass ratio 3:1 including both ini- 
tially nonspinning BHs (the "A" cases) and BHs with 
spin parameter set to 0.75 initially (the "B" cases). Ta- 
ble [III specifies the AMR grid structure used in the simula- 
tions and Table HITl summarizes the quantitative results. 
For readers interested only in a brief summary of the 
most interesting results, please skip to Sec. |Vl Detailed 
simulation results are described below. 



A. Magnetic Field Study: Nonspinning Black Hole 

Figure [2] shows density contours on the orbital plane at 
selected times for the unmagnetized, zero BH spin case, 
AO. Notice that the NS density contours in the top- left 
plot are nearly unchanged after three orbits (top-center 
plot), confirming that the initial data are consistent with 
quasiequilibrium. After about 3.5 orbits the NS tidally 
disrupts (top-right plot), and about 95% of the NS mat- 
ter promptly falls into the BH. Matter in the low-density 



NS outer layers far from the accreting funnel of mat- 
ter forms a tidal tail that wraps around the black hole 
and smashes into itself near the BH, generating a large 
amount of shock heating (bottom- left plot). Meanwhile, 
accretion slows considerably. After intersecting itself, the 
inner regions of the tidal tail forms a disk that orbits 
the BH, while the fluid velocity distribution in the outer 
tail indicates slow accretion onto the disk (bottom- middle 
plot). Shortly after disk formation, only about 2% of the 
NS matter remains outside the BH (bottom-right plot), 
and the density directly outside the AH begins to plum- 
met, ultimately forming a low-density cavity around the 
BH similar to the one shown more prominently in the 
bottom left frame of Fig. [5l This cavity indicates the 
presence of an innermost stable circular orbit (ISCO). 

Figure [3] shows the accretion history for all cases in 
which the BH has zero spin initially. Regardless of mag- 
netic field configuration or strength, hy t ^ 570M about 
95% of the NS matter - including the most strongly mag- 
netized matter in the star - has been accreted by the BH. 
After this violent merger, the only case that noticeably 
deviates from the magnetic-free case (AO) is case A4, the 
case in which the initial seed magnetic fields are both 
strong (|5|max ^ lO^^G) and pushed to the NS surface 
(Pc = 0.001). The disk mass in case A4 is two times 
larger than any other case, but the final disk is only about 
2.5% of the initial NS rest mass - much smaller than in 
similar BHNS simulations with a moderate aligned BH 
spin, in which disk masses of ~ 10% are common. Case 
A3 is identical to case A4, except for the fact that its 
seed magnetic fields are about an order of magnitude 
weaker (|P|max ^ lO^^G), and its accretion history is 
virtually indistinguishable from that of case AO. Thus it 
is the magnetic field strength - and not the different ge- 
ometries explored here - that significantly infiuences the 
dynamics. The final accretion rate implies a disk half-life 
of between 3500-5500M or lOO-15O(Mo/1.4M0)ms (de- 
pending on what points are chosen to calculate the final 
slope). 

Figure 2] plots magnetic energy outside the AH versus 
time, for all magnetized nonspinning BH cases studied. 
Magnetic fields were added shortly before tidal disrup- 
tion. The magnetic fields do not change significantly in 
the NS prior to disruption, but at the point of disruption, 
there are two competing effects that infiuence the mag- 
netic energy. For one, the NS is being tidally disrupted, 
stretching the magnetic field lines, amplifying the mag- 
netic field strength and energy. On the other hand, the 
magnetized fiuids comprising the NS are being rapidly 
accreted into the BH. Our results show that there is a 
slight amplification of magnetic energy during tidal dis- 
ruption, but after tidal disruption the magnetic energy is 
always less than the magnetic energy in the seed magnetic 
fields. The evolution is followed for ~ 5O(Mo/1.4M0)ms 
(assuming NS rest mass of I.4M0) after disk formation, 
but ultimately no large magnetic energy amplification is 
observed. This is likely due to the fact that the magnetic 
fields in the disk are mostly toroidal (see Fig. [5|) and once 
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TABLE II. Grid configurations. Here, A^ah denotes the number of grid points covering tfie diameter of the (spherical) AH 
initially, and A^ns denotes the number of grid points covering the smallest diameter of the neutron star initially. 

Case Grid Hierarchy (in units of M)*^"-* Max. resolution A^ah A^ns 

A0-A4 (196.7, 98.35, 49.18, 24.59, 12.29, 6.147, 3.073, 1.414 [1.660]) M/32.5 41 85 

B0-B4 (210.2, 92.49, 46.24, 23.12, 11.56, 5.780, 2.890, 1.445 [1.642], 0.7554 [N/A]) M/60.9 56 80 

^"'^ There are two sets of nested refinement boxes: one centered on the NS and one on the BH. This column specifies the half 
side length of the refinement boxes centered on both the BH and NS. When the side length around the NS is different, we 
specify the NS half side length in square brackets. If there is no corresponding NS refinement box (as is the case when the NS 
is significantly larger than the BH), we write [N/A] for that box. 



TABLE HI. Magnetized BHNS simulation results. Here Mdisk is the rest mass of the material outside the AH at the end of the 
simulation, af = clbu/Mbu is the spin parameter of the BH at late times determined by the isolated horizon formalism. The 
total energy and angular momentum carried off by the gravitational radiation are given by AEqw and AJgw, respectively. 
Vkick is the kick velocity due to recoil. A^orbits specifies the number of orbits in the inspiral phase before merger, defined as 
the time at which the (2,2) mode of the GW amplitude reaches maximum. SE and SJ measure the violation of energy and 
angular momentum conservation, as defined in Eqs. (|26p and (|27p . respectively, at the end of the simulation. The error in Vkick 
is estimated by comparing the results obtained by several GW extraction radii between 50M-100M. 



Case 


Mdisk/Mo 


af 


AEgw/M AJgw/J 


^;kick (km/s) 


A^orbits 


6E 


6J 


AO 


0.019 


0.55 


0.011 


0.20 


40 ±2 


4.8 


0.2% 


2% 


Al 


0.017 


0.55 


0.011 


0.20 


40 ±2 


4.8 


0.2% 


2% 


A2 


0.017 


0.55 


0.011 


0.20 


40 ±2 


4.8 


0.2% 


2% 


A3 


0.015 


0.55 


0.011 


0.20 


40 ±2 


4.8 


0.2% 


3% 


A4 


0.028 


0.55 


0.011 


0.20 


40 ±2 


4.8 


0.1% 


1% 


BO 


0.098 


0.84 


0.011 


0.15 


67 ±6 


6.9 


0.6% 


8% 


Bla 


0.090 


0.84 


0.011 


0.15 


67 ±6 


6.9 


0.6% 


8% 


Bib 


0.090 


0.84 


0.011 


0.15 


67 ±6 


6.9 


0.6% 


8% 


B2 


0.090 


0.84 


0.011 


0.15 


67 ±6 


6.9 


0.6% 


7% 


B3 


0.090 


0.84 


0.011 


0.15 


67 ±6 


6.9 


0.6% 


7% 


B4 


0.117 


0.85 


0.010 


0.14 


54 ±4 


6.8 


0.4% 


7% 



the disks have formed magnetic winding saturates. Am- 
plification of magnetic fields by instabilities such as the 
magnetorational instability (MRI) may occur, but the 
resolution in our simulations may not be high enough to 
resolve the small-scale turbulence associated with these 
instabilities. 

Notice that the magnetic energy in case A2 is about 
three times that of Al, both initially and when the simu- 
lation was terminated. These cases differ only in the de- 
gree of central condensation of the initial magnetic field 
(see Table U). 

Cases A2 and A3 are identical except for the pressure 
cutoff of the seed magnetic fields. In case A2 magnetic 
fields are set to zero for pressures P < Pc = 0.04Pmax, 
where Pmax is the maximum pressure of the NS. However, 
in case A3 Pc is set to 0.001, so the seed magnetic fields 
are pushed much closer to the NS surface. This results in 
about a 16% amplification of initial magnetic energy in 
case A3 (Fig. |4]), but increases the final magnetic energy 
by about a factor of 15. This is consistent with the fact 
that the core of the NS is invariably accreted into the BH 
during merger in these BHNS simulations, so the disk is 
comprised of what were the outer layers of the NS. Thus, 



the stronger the seed magnetic field in the outer layers of 
the NS, the stronger and more dynamically relevant the 
magnetic fields in the disk. 

Cases A3 and A4 differ only in initial seed magnetic 
field strength; the seed magnetic fields are uniformly 
about an order of magnitude stronger in case A4. Fig- 
ure [3] demonstrates that the physical extent of the disk 
is very strongly influenced by the strong seed magnetic 
fields of case A4. Figure |4] reinforces that observation; 
though only about 1% of the NS rest mass exists in the 
disk of A3, less than 0.5% of the seed magnetic energy 
remains in the disk. Compare this to case A4, where the 
disk mass is about twice as large, but where more than 
an order of magnitude more magnetic energy remains in 
the disk. 

Magnetic fields play an important dynamical role in 
only one nonspinning case, A4, amplifying the disk mass 
by a factor of two. Figure [5] shows how the magnetic 
field configuration evolves in this case, from t = 448. 5M 
when the magnetic fields are first seeded into the NS 
(top-right), until the simulation is stopped at t = 2620M 
(bottom-right). Magnetic field lines are greatly stretched 
during disk formation (third row on the right), result- 
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FIG. 2. Orbital-plane rest mass density contours at selected times for case AO. Contours are plotted according to po — 
Po,max(10~°"^^^), (j=0, 1, ... 5), with darker greyscaling for higher density. The maximum initial NS density is hvpo,ma^ — 0.126, 
or y9o,max = 9 X lO^^g cm~^ (IAMq / Mq^ . Arrows represent the velocity field in the orbital plane. The black hole AH interior 
is marked by a filled black circle. The ADM mass for this case is M = 2.5 x lO~^(Mo/1.4M0)s= 7.6(Mo/1.4M0)km. 



ing in a strongly-magnetized disk. At late times, the 
magnetic fields that remain in the disk are very tightly 
wound. The bottom left and center frames of the figure 
show that a cavity has formed around the BH near the 
end of the simulation. This hollow region is indicative of 
the presence of an ISCO, which is interesting because it 
has been suggested that stresses in magnetized disks may 
suppress the presence of an ISCO [96|, |97[. However, it 
may be that longer simulations are required for the disk 
cavity to be filled. 

The significant boost in disk mass in case A4 indicates 
that with sufficiently strong NS magnetic fields, merger 
dynamics may be significantly affected. Figure [6] com- 
pares orbital-plane density contours in cases AO (unmag- 
netized) and A4 during the late stages of tidal disrup- 
tion. Notice that the stronger magnetic pressures in case 
A4 push out the outer layers of the NS during tidal dis- 
ruption. Nevertheless, Fig. [7] shows that the magnetic 
pressure in this matter distribution does not exceed 3% 
the gas pressure at the same time as Fig. [6l Therefore, 
large changes in remnant disk mass do not require huge 
magnetic-to-gas pressure ratios. 

Figure [8] compares AO and A4 rest-mass density pro- 



files when these simulations were stopped at t = 2620M. 
The density profiles in these two cases are similar; at the 
end of the simulation, low-density matter still flows into 
the BH from the poles. For A4, it is clear from the 6^ 
contour plot (bottom graph) that most of the magnetic 
fields are confined inside the remnant disk near the equa- 
torial plane, consistent with the magnetic field-line plots 
in Fig. El 

Compared to case AO, the final disk in case A4 is more 
massive, and the accretion rate is lower (Fig. [3]). Fig- 
ure [9] compares the distribution of rest-mass density and 
entropy (log K) in the remnant disks of cases AO and A4, 
at the time shortly after disk formation. The disk in case 
AO is of roughly uniform rest-mass density in the orbital 
plane. Though the disk volume in cases A4 and AO are 
comparable, the case A4 disk is about twice as massive 
as AO. Entropy in case A4 is lower, and entropy contours 
fall off rapidly with density. On the other hand, the en- 
tropy in case AO is much more uniform in the disk, with 
only a slight drop near the BH. Thus in the nonspinning 
case, adding strong seed magnetic fields to the NS results 
in colder, denser, more massive disks. 

Figure [10] shows I? profiles for case A4's disk, plot- 
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FIG. 3. Rest-mass fraction outside the AH for all cases in 
which the BH is initially nonspinning. 



Short-wavelength variations in 1? contours appear near 
the BH only, with longer-wavelength variations outside. 
This is likely a numerical artifact, since the disk spans 
about four AMR refinement levels, which are centered 
on the BH. Each refinement level drops the resolution by 
a factor of two, meaning that matter in the outer reaches 
of this disk is more poorly resolved by a factor of 16 than 
the region near the BH. This filtering of wavelengths due 
to AMR likely suppresses magnetic-induced turbulence 
in the disk. 

Figure [TTl contrasts 6^/(2P) contours for cases A3 and 
A4. As 6^/(2P) approaches unity, magnetic fields should 
have more of a dynamical impact. Case A3 is identical 
to A4, except its seed magnetic fields are weaker by an 
order of magnitude. Correspondingly, typical values of 
Ip' / (2P) in this late-time disk are roughly 2-3 orders of 
magnitude lower in case A3 than A4. This is consistent 
with the finding that the remnant disk mass and accre- 
tion rate are unaffected by magnetic fields in case A3; 
pressure from the magnetic fields is simply too small to 
be dynamically relevant. Thus, initial internal magnetic 
pressures of order 0.1% of gas pressure may be required 
for magnetic fields to have a significant impact on the 
system's dynamics. 
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FIG. 4. Total magnetic energy M outside the AH for all 
magnetized cases in which the initial BH possesses zero spin, 
normalized by the ADM mass M = 9.3x lO^^(Mo/1.4M0)erg. 



ted at the same time as in Fig. [9l The strong magnetic 
fields in A4 enhance the final disk mass by more than a 
factor of two, even though magnetic pressure never ex- 
ceeds about 3% of the gas pressure, with typical values 
around 0.01% (left plot)-about an order of magnitude 
lower than during tidal disruption (Fig. [71 left frame). 



B. Magnetic Field Study: Spinning Black Hole 

Figure [T2l outlines the basic evolution scenario for the 
fiducial aBH/^BH=0.75 (henceforth "spinning") case, 
BO. Unlike the nonspinning case, the NS is slightly dis- 
torted from its equilibrium shape after about three or- 
bits (top- left plots, cf. top- left plots in Fig. [2]). Although 
the spinning cases start with the same initial orbital an- 
gular frequency as the nonspinning cases, all spinning 
cases require about two more orbits before an accretion 
funnel forms (top-right frame). This is due to the well- 
known "orbital hang- up" effect. Unlike the nonspinning 
case, in which about 95% of the NS matter immediately 
funnels into the BH, strong frame-dragging in the spin- 
ning cases twists the accreting funnel around the BH, 
promptly accreting only ~70% of the NS matter. Ac- 
cretion slows after the funnel twists around the BH and 
intersects itself, generating shock heating and produc- 
ing a small, dense disk-like structure that expands and 
rarefies as it orbits the BH. Notice that this "disk" is 
much smaller and denser than in the nonspinning cases 
(bottom- left plot, cf. Fig. [2]). Attached to this small 
"disk" , the outer layers of the NS have formed a long tidal 
tail, which is ejected to a large radius before slowly falling 
back onto the expanding but accreting "disk" (bottom- 
middle frame). The disk continues its expansion as much 
of the tail falls into it. Near the time when the simu- 
lation is stopped, we find no indication of a cavity near 
the BH (bottom-right frame). Longer simulations may 
be required for a cavity to appear, which might indicate 
the presence of an ISCO. Note that in contrast with the 
expected long-term evolution of this unmagnetized case. 
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FIG. 5. Orbital-plane rest-mass density contours (left column), 3D density profiles (middle column) and 3D magnetic field lines 
(right column) at four selected times for case A4. The times in the rows (top to bottom) are t/M =448.5, 515.4, 613.7, and 2620. 
Density contours in the orbital plane (left column) are plotted according to po = po,max(10~°'^^-^), (i=0, 1, ... 5), with darker 
greyscaling for higher density. The maximum initial NS density is Atpcmax = 0.126, or /?o,max = 9 x lO^^g cm~^(1.4M0/Mo)^. 
Arrows in density contour plots represent the velocity field in the orbital plane, and the black hole AH interior is marked by a 
filled black circle. Magnetic fields are plotted as streamlines of the magnetic field vector B*, distributed in proportion to \B'^\. 
The ADM mass for this case is M = 2.5 X lO~^(Mo/1.4M0)s= 7.6(Mo/1.4M0)km. The 3D visualizations were produced using 
the ZIBamira software system [oBI ]. 
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FIG. 6. Rest-mass density and velocity profile snapshots during NS tidal disruption for cases AO (left) and A4 (right). Density 
contours are plotted in the orbital plane according to po = po,max(10~°"^^'^ ), (i=0, 1, ... 5), with darker greyscaling for higher 
density. The maximum initial NS density is npo^ma^ = 0.126, or po,max = 9 x lO^^g cm~'^(1.4M0/Mo)^. Arrows represent the 
velocity field in the orbital plane, and the black hole AH interior is marked by a filled black circle. The ADM mass for this 
case is M = 2.5 x lO"^(Mo/1.4M0) s= 7.6(Mo/1.4M0)km. 




FIG. 7. Pressure ratio 6^/(2P) (left) and magnetic pressure 6^/2 (right) contours during NS tidal disruption, at the same 
time as Fig. [1 plotted according to b'^/{2P) = 10"^-^(10"^-^^), (j=0, 1, ... 5), and nb'^ = 10"^(10"^-^^ ), (j=0, 1, ... 5). 
Darker greyscaling denotes higher values. Contours are only plotted for regions with densities higher than the lowest- density 
po contours in Fig. H In cgs units, ^ 6 x lO^^dyn cm~^(1.4M0/Mo)^. 
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FIG. 8. Rest-mass density contours for case AO (top), A4 
(middle), and magnetic energy density (bottom) profile 
for A4 in the meridional plane at the end of simulations 
(t = 2620M) . Density contours are plotted according to po = 
yOo,maxlO~^'^^°'^-^^-^ , (j = 0, ...,6), with darker greyscaling for 
higher density, contours are plotted according to nb^ — 
^Q-i2+o.833j^ (j = 0, 6). The maximum initial NS density 

is /€po,max = 0.126, or po,max = 9 X lO^^g ciR-^ {1 AMq / Mof . 

Arrows represent the velocity field in the meridional plane. 
In cgs units, the ADM mass for this case is M = 2.5 x 
lO-^(Mo/1.4M0) s= 7.6(Mo/1.4M0)km, and k'^ = 6 x 
lO^^erg cm-^(1.4M0/Mo)^ = 7 x lO^^g cm-2(1.4M0/Mo)^ 



stresses in magnetized disks may suppress the formation 
of hollow regions around the ISCO [96, 97]. 

The corresponding accretion history for all cases in 
which the BH initially has spin is shown in Fig. [131 Simi- 
lar to the nonspinning cases, only the most strongly mag- 
netized case, B4, has an appreciably different accretion 
history. The NS in case B4 possesses the same magnetic 
field geometry and strength as the NS in the (nonspin- 
ning) case A4: the seed magnetic fields are both strong 
(1^1 max ^ lO^^G, initially) and pushed all the way to 
the NS surface (Pc = 0.001). The magnetic fields in 



case B4 increase the disk mass from about 10% to 14%, 
a net 40% amplification in disk mass. Notice also that 
BH spin alone has a very significant influence on final 
disk mass; increasing initial aligned BH spin from zero 
to aBH/^BH=0.75 increases the final disk mass by about 
an order of magnitude (~0.9% - '^10%). Based on the 
final accretion rate, the half-life of the disk is roughly 
5000M , or - 14O(Mo/1.4M0)ms. 

Figure [TH plots the magnetic energy outside the AH 
versus time for all spinning, magnetized cases. During 
merger there are two competing effects: magnetic energy 
will increase as the NS tidally disrupts and the field lines 
are stretched, while magnetic energy outside the AH will 
decrease as magnetized NS matter is accreted into the 
BH. This major accretion event occurs at t ~ 900M for 
all spinning cases, corresponding to a spike in magnetic 
energy at that time. At t = lOOOM, only about 20% of 
the NS matter remains outside the BH, but the magnetic 
energy in all cases is higher than when the seed fields were 
added to the NS, indicating that a large amplification in 
magnetic field strength has occurred. Only case B4 has 
sufficiently strong magnetic fields to severely impact the 
disk dynamics; the disk mass is amplified by about 40%, 
and the magnetic energy in this case is 1-2 orders of 
magnitude stronger than any other case at all times. 

To determine the sensitivity of our results on the time 
at which the magnetic fields were added to the NS, two 
simulations were performed: cases Bla and Bib. These 
simulations are identical except the seed magnetic fields 
were added to the NS about half an orbit later in case 
Bib. The accretion histories of cases Bla and Bib over- 
lap completely, implying that the bulk dynamics are un- 
affected by when the magnetic fields were inserted into 
the NS. During and directly after the merger {t ~ 900M), 
magnetic energy in cases Bla and Bib overlap. After 
merger, only a tiny fraction of the fields in the outer lay- 
ers of the NS remains outside the AH. Correspondingly, 
the magnetic energy plummets in both cases to negligibly 
small values, yet agree to within an order of magnitude 
even at late times. Thus we conclude that the final re- 
sult is largely insensitive to when the seed magnetic fields 
were added. 

Next we analyze how magnetic fields evolve over time 
in case B4, starting with the time at which they were 
seeded into the NS (top-right frame of Fig. [T5]) . At the 
onset of tidal disruption (second row), the magnetic field 
structure has changed significantly, even within the non- 
disrupted regions of the NS. Apparently the magnetic 
fields have undergone some slight rearrangement since 
they were added to the NS. 

After the accretion funnel has wrapped around the BH 
and intersected itself, it forms a small disk-like structure 
around the BH (bottom frames). In this "disk" region, 
the magnetic fields wind around the BH. The B4 sim- 
ulation is continued for about 3O(Mo/1.4M0)ms after 
tidal disruption, and then it is stopped. At this time, 
the BH+disk system has drifted significantly (bottom 
plots of[T5|), due to gauge effects and the gravitational- 



14 




-15 -10 -5 5 10 15 -15 -10 -5 5 10 15 

X / M X / M 

FIG. 9. Top two plots: Rest-mass density and velocity snapshots shortly after disk formation for cases AO (left) and A4 (right). 
Density contours are plotted in the orbital plane according to po = yOo,max(10~°"^^"' ), (j—0, 1, ... 5), with darker greyscaling 
for higher density. The maximum initial NS density is Atpcmax = 0.126, or po,max = 9 x lO^^g cm~^(1.4M0/Mo)^. Arrows 
represent the velocity field in the orbital plane. The black hole AH interior is marked by a filled black circle. The ADM mass 
for this case is M = 2.5 x lO"^(Mo/1.4M0) s= 7.6(Mo/1.4M0)km. 

Bottom two plots: Snapshots of the entropy parameter K contours for cases AO (left) and A4 (right). The light grey regions 
correspond to 1.4 < log^^Q K < 2.6, and the dark grey region corresponds to 2.6 < log^^Q K < 3.8. 



wave kick at merger. There is a great deal of winding 
of magnetic fields threading the disk at this time (green 
lines), but no strong evidence of collimation around the 
BH poles or magnetic field turbulence in the disk. The 
lack of magnetic field turbulence may be due to in- 



sufficient resolution in the disk, which artificially sup- 
presses instabilities like MRL Insufficient resolution, cou- 
pled with the termination of the simulation after only 
3O(Mo/1.4M0)ms may explain why no magnetic field 
collimation was observed. The bottom left and center 
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FIG. 10. h^/(2P) (left) and I? (right) contours, at the same time as Fig. [1 plotted according to h^/(2P) = 10~^-^(10"^-^^ ), 
(j=0, 1, ... 5), and Kh^ = 10"^(10"^-^^), (j=0, 1, ... 5). h^/{2P) and contours are only plotted for regions with densities 
higher than the lowest- density po contours in Fig. (9] In cgs units, hi~^ = 6 x lO'^^erg cm~'^(1.4M0/Mo)^. 




FIG. 11. h^/{2P) contours for cases A3 (left) and A4 (right), plotted according to h^/{2P) = 10"^-^(10"^-^^), (j=0, 1, ... 5). 
Contours are only plotted for regions with densities above the low-density cutoff in Fig. [9l 



frames of the figure show that a cavity has not formed 
around the BH by the end of the simulation. This result 
appears to be consistent with studies which suggest that 
stresses in magnetized disks may suppress the presence 
of an IS CO [9fi,^]. However, longer, more accurate disk 



evolutions will be necessary to fully assess the agreement 
between these studies and simulations in full GR. 

The strong seed magnetic fields in case B4 amplify the 
disk mass significantly, similar to case A4. Figure [16] 
demonstrates how the NS density contours are affected 
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FIG. 12. Orbital-plane density contours at selected times for case BO. Density contours are plotted in the orbital plane, 
according to po = po,max(10~°"^^-^ ), (j=0, 1, ... 5), with darker greyscaling for higher density. The maximum initial NS 
density IS /x^PCmax = 0.126, or /9o,max — 9 X lO^^g cm ^ [IAMq /Mq)'^ . Arrows in density contour plots represent the velocity 
field in the orbital plane, and the black hole AH interior is marked by a filled black circle. The ADM mass for this case is 
M = 2.5 X lO"^(Mo/1.4M0)s= 7.6(Mo/1.4M0)km. 



by such strong magnetic fields, comparing case B4 (right 
plot) with the unmagnetized case (left plot) shortly af- 
ter tidal disruption. Although the seed magnetic fields 
are the same strength and geometry as in case A4 (cf. 
Fig. [6|), the NS outer layers in case B4 are pushed out 
much more than in case A4. The magnetic pressure in 
case B4 does not rise above about 3% of the gas pres- 
sure in the orbital plane, as shown in Fig. [T71 (left plot). 
This is consistent with case A4. Although 6^/ (2P) varies 
strongly throughout the NS, 6^ tends to be stronger in 
the densest regions of the NS (right plot), again similar 
to case A4 (cf. Fig. fTO]) . 

Figure [T8l contrasts the distribution of rest-mass den- 
sity and entropy (logK) in the remnant disks of cases 
BO and B4, at the time in which the B4 simulation was 
stopped. In the nonspinning cases AO and A4 (Fig. [9]), 
the sizes of the remnant disks in the orbital plane are 
remarkably similar, though the disk mass in case A4 is 
about twice AO's disk mass. The extra disk mass in case 
A4 may be explained by the existence of a high-density 
ring of matter in the disk. Unlike the nonspinning cases 
however, the high-density regions of the remnant disks 



in cases BO and B4 are remarkably similar (upper plots). 
Despite this, B4's final disk possesses about 40% more 
mass than BO. Apparently the excess mass in case B4's 
disk comes from its larger volume. The size difference be- 
tween BO and B4's final disks is seen more clearly in the 
entropy contour plots (bottom two plots), which demon- 
strate the disks are hotter close to the BH and in the 
lowest-density outer regions. 

Magnetic field amplitude contour plots for case B4's 
disk at the same time as Fig. [18] are displayed in Fig. [T9l 
As in case A4, the final disk mass is greatly increased by 
strong magnetic fields, even though magnetic pressure 
never exceeds about 3% of the gas pressure, with typ- 
ical values around 0.1% (left plot of Fig. [T9|) . Bubbles 
of enhanced 6^, corresponding to order-of- magnitude in- 
creases in appear in a small ring around the BH and 
in the disk's outer regions (right plot of Fig. [T9|) . 

Finally, Fig. [20] compares BO and B4 profiles in the 
meridional plane. As in the nonspinning cases, the ge- 
ometry of the disk appears to be mostly unchanged by 
the addition of magnetic fields. There is small amount 
of material fiowing into the BH from the poles, and the 



17 



10-3 ^ ^ ^ ^ ^ ^ ^ ^ , 




0.05 ' i i i I i i i I I I I I I I I 

500 1000 1500 2000 

t/M 

FIG. 13. Rest-mass fraction outside the AH for all 
cases in which the initial BH possesses spin parameter of 
aBH/MBH=0.75. 



magnetic fields are mainly confined inside the disk. 



C. Gravitational Waves Study 

Gravitational waves are extracted at several radii be- 
tween 50M and 130M. These radii are sufficiently far 
from the binary for the waves to overlap very well when 
plotted against retarded time, after accounting for the 
amplitude fall-off with radius. The upper panel of Fig. [211 
shows the dominant (2,2) mode of the gravitational wave- 
form as a function of retarded time tret for case AO, and 
the difference between the two data sets is displayed in 
the lower panel. The difference in amplitude between 
AO and A4 waveforms is less than 3%, indicating that 
magnetic fields explored here do not significantly impact 
the global dynamics during BHNS inspiral and merger. 
By contrast, the difference between the BO and B4 wave- 
forms during merger is substantial (Fig. [22]) . The more 
highly spinning BHs of cases BO and B4 possess smaller 
ISCOs, enabling NS to orbit the BH more closely before 
accreting, resulting in more gravitational wave cycles. 
Further, at smaller separations from the BH, the effects 
of frame dragging are much more pronounced, twisting 
NS matter around the BH and reducing the BH accretion 
rate, ultimately giving the magnetic fields more time to 
amplify and affect the global dynamics of NS matter. 

However, the observability of magnetic effects on grav- 
itational waveforms depends on the response of the GW 
detectors to the time series data. One way to assess the 
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FIG. 14. Total magnetic energy outside the AH for all 
magnetized cases in which the initial BH possesses spin pa- 
rameter of aBH/MBH=0-75, normalized by the ADM mass 
M = 9.3 X lO^^(Mo/1.4M0)erg.. 

detectability is to compute the mismatch 

MM ^ 1 - jhsolhs,) _ ^28) 

V{hBo\hBo){hB4\hBA) 

between the waveforms of BO and B4, assuming the NS 
mass is 1.4 Mq, where 

Here h = h-^ — ihx, h is the Fourier transform of h{t), 
and Shif) is the power spectral density of the Advanced 
LIGO noise. Using the Advanced LIGO broadband con- 
figuration HIGH_DET_HIGH_P, the minimum mismatch 
(by varying the time and phase shifts between the two 
waveforms) between BO and B4 waveforms is only 0.004, 
indicating that it may be challenging for Advanced LIGO 
broadband to detect the strong internal magnetic fields 
of case B4. 

Complete GW spectra in the frequency domain require 
the creation of hybrid waveforms, which stitch together 
numerical and post-Newtonian (PN) waveforms. We gen- 
erate hybrids by first computing the minimum of 

pdt I (hf^-^rr (30) 

]l ie{+,x} 

via the Nelder-Mead algorithm, using as free parame- 
ters initial PN phase, amplitude, and orbital angular 
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FIG. 15. Orbital-plane density contours (left column), 3D density profiles (middle column) and 3D magnetic field lines (right 
column) at four selected times for case B4. The times in the rows (top to bottom) are t/M =634.8, 744.1, 931.2, and 2003. 
Density contours in the orbital plane (left column) are plotted according to po = po,max(10~°'^^-^), (i=0, 1, ... 5), with darker 
greyscaling for higher density. The maximum initial NS density is Atpcmax = 0.126, or /?o,max = 9 x lO^^g cm~^(1.4M0/Mo)^. 
Arrows in density contour plots represent the velocity field in the orbital plane, and the black hole AH interior is marked by a 
filled black circle. Magnetic fields are plotted as streamlines of the magnetic field vector B*, distributed in proportion to \B'^\. 
The ADM mass for this case is M = 2.5 X lO~^(Mo/1.4M0)s= 7.6(Mo/1.4M0)km. The 3D visualizations were produced using 
the ZIBamira software system [95l ]. 
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FIG. 16. Density and velocity profile snapshots during NS tidal disruption for cases BO (left) and B4 (right). Density contours 
are plotted in the orbital plane according to po = po,max(10~°"^^-^), (i=0, 1, ... 5), with darker greyscaling for higher density. 
The maximum initial NS density = 0.126, or po,max — 9 X lO^'^^g cm ^ (IAMq /Mq)'^ . Arrows represent the velocity 

field in the orbital plane, and the black hole AH interior is marked by a filled black circle. The ADM mass for this case is 
M = 2.5 X lO"^(Mo/1.4M0) s= 7.6(Mo/1.4M0)km. 




FIG. 17. Pressure ratio /(2P) (left) and magnetic pressure 6^/2 (right) contours during NS tidal disruption, at the same 
time as Fig. [H plotted according to b'^/{2P) = 10"^-^(10"^-^^'), (j=0, 1, ... 5), and nb'^ = 10"^(10"^-^^'), (j=0, 1, ... 5). 
Darker greyscaling denotes higher values. Contours are only plotted for regions with densities higher than the lowest- density 
po contours in Fig. [H In cgs units, ^ 6 x lO^^dyncm~^(1.4M0/Mo)^. 
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FIG. 18. Top two plots: Density and velocity profile snapshots at the time in which the BO simulation is stopped, for cases BO 
(left) and B4 (right). The contours represent the density in the orbital plane, plotted according to po = yOo,max(10~°"^^-^), 
(j=0, 1, ... 5), with darker greyscaling for higher density. The maximum initial NS density is K,po,max = 0.126, or po,max = 
9 X lO^^g cm~^(1.4M0/Mo)^. Arrows represent the velocity field in the orbital plane. The black hole AH interior is marked 
by a filled black circle. The ADM mass for this case is M = 2.5 x lO~^(Mo/1.4M0)s= 7.6(Mo/1.4M0)km. 
Bottom two plots: Snapshots of entropy parameter K contours for cases BO (left) and B4 (right). The light grey regions 
correspond to 1.4 < log^^Q K < 2.6, and the dark grey region corresponds to 2.6 < log^^Q K < 3.8. 



frequency. Here, and h^^^ specify our numeri- 

cal BH NS waveforms and the TaylorTl PN waveforms 
of |lOO| . respectively. The integration bounds were cho- 
sen to be ti ^ 200M and U ~ 400M. The hybrid wave- 
form consists of a linear combination of the PN and NR 



waveforms, as in [101|. 

Effective GW strains of the hybrid waveforms in fre- 
quency space are plotted in Fig. [23] for AO, and Fig. [23] 
for cases BO and B4. Assuming the NS has a rest 
mass of 1.4M0 and binary distance of lOOMpc, we plot 




FIG. 19. Pressure ratio b^/{2P) (left) and magnetic pressure 6^/2 (right) contours, at the same time as Fig. [181 plotted 
according to b'^/{2P) = 10"^-^(10"^-^^), (j=0, 1, ... 5), and nb'^ - 10"^(10"^-^0, (j=0, 1, ... 5). b^/{2P) and 6^ contours 
are only plotted for regions with densities higher than the lowest- density po contours in Fig. 1181 In cgs units, Hi~^ = 6 x 
lO^^dyn cm-2(1.4M0/Mo)^ 



these strains against the Advanced LIGO sensitivity 
curve /iLiGo(/) = \/ fSh{f)' Within this distance, the 
BHNS event rate is estimated to be 5 x 10~^-10 per 
year, assuming an overall rate of 0.05-100 mergers per 
Myr per Milky Way-equivalent galaxy (and a density of 
0.1 gal/Mpc^) [29|. Cases A1-A4 and B1-B3 are not 
shown in the figures because their effective GW strains 
are almost indistiguishable from cases AO and BO, re- 
spectively. Advanced LIGO in the chosen configuration 
may be able to marginally distinguish between BHNS 
waveforms and those produced by BHBH mergers at high 
frequencies (500-lOOOHz). However, the effects of even 
the strongest magnetic fields chosen here (cases A4 and 
B4) on the waveforms are quite small and may be chal- 
lenging for Advanced LIGO to detect in a broadband 
configuration. This is consistent with the result from 
the minimum mismatch calculation between BO and B4, 
as mentioned above. Nevertheless, recent innovations ex- 
ploring 'squeezed light' effects may reduce quant um n oise 
and increase the sensitivity in this very domain |l02| . 

It has been suggested that the differences between 
BHNS and BHBH waveforms during late inspiral and 
merger may be used to extract the tidal deformability 
of the NS [103], which can in turn be used to constrain 
the NS EOS. Our results suggest that seeding the NS 
with magnetic fields of the configurations and strengths 
explored here do not alter gravitational waveforms sig- 
nificantly. Further studies at higher resolution may be 
required to confirm this finding. 



D. Energy and angular momentum conservation 

We compute the energy AEqw and angular momen- 
tum A Jgw carried away by the GWs, as weh as the GW 
recoil velocity I'kick- Violation of energy 5E and angular 
momentum (5 J, as defined in Eqs. (|26|) and (|27|) , respec- 
tively, is also monitored. Results are given in Table lllli 
For cases A0-A4, SE ^ 0.2% and SJ ^ 2%. Whereas for 
cases B0-B4, SE ^ 0.6% and 5 J ^ 8%. In all cases, a 
fraction of E and J are lost spuriously, and the situation 
is slightly worse in spinning BH cases. 

To further study the issue of E and J loss, we calculate 
the angular momentum inside the computational domain, 
Jint, using Eq. ([22]) . as well as the accumulated angular 
momentum carried away by the GWs, A Jew- Figures [251 
and [26l show the evolution of Jint and AJqw for represen- 
tative cases (A4 and B4). The corresponding plots for the 
various components of energy are similar. Conservation 
of angular momentum implies that Jsum = Jint A Jew is 
constant in time and is equivalent to the ADM angular 
momentum of the binary. Notice that the total angu- 
lar momentum Jsum is conserved well during inspiral for 
both cases (tret < 500M for case A4 and tret < 700M for 
B4). Post-merger, A Jew flattens as GW emission sub- 
sides in all cases. Jint also flattens after merger in case 
A4, as expected. Hence the spurious loss of J 2%) 
occurs primarily during merger for case A4. The same 
behavior is observed in cases AO- A3 as well. However, in 
cases B0-B4, Jint decreases on a secular timescale after 
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FIG. 20. Rest-mass density profile for case BO (top), B4 
(middle) and magnetic energy density (bottom) profiles 
for B4 in the meridional plane at the end of simulation (t = 
2003. 2M). Density contours are plotted according to po = 
po,maxlO~^"^^°"^"^^^ , (j = 0, ...,6), with darker greyscaling for 
higher density, contours are plotted according to nb^ = 
^Q-i2+o.833j^ (j = 0, ...,6). The maximum initial NS density 
is Atpcmax = 0.126, or po,max = 9 X lO^^g cui-^ {1 AMq / Mq)'^ . 
Arrows represent the velocity field in the meridional plane. 
In cgs units, the ADM mass for this case is M = 2.5 x 
lO-^(Mo/1.4M0) s= 7.6(Mo/1.4M0)km, and k'^ = 6 x 
lO^^erg cm-^(1.4M0/Mo)^ = 7 x lO^^g cm-2(1.4M0/Mo)^ 
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FIG. 21. Upper panel: (2,2) mode of the gravitational wave 
strain (solid line) and /i22 (dashed line) versus retarded 
time tret for case AO. 

Lower panel: difference in (solid line) and /i22 (dashed 
line) between cases AO and A4 versus tret , normalized by the 
maximum value of \h22 — '^^22! ^^r case AO over time, /i22max. 



lution. 

Figures [25l and [26l also show the evolution of the BH's 
angular momentum Jbh, computed using the isolated 
and dynamical horizon formalism. Jbh monotonically 
increases in both A4 and B4, and is always less than 
Jint- This seems to suggest that the spurious loss of J 
in case B4 after merger may occur mainly in the rem- 
nant disk, where a substantial amount of matter is con- 
tinuously crossing the refinement levels. However, the 
possibility of spurious loss of J in the strong-field region 
near the BH cannot be ruled out completely. One way 
to resolve the issue would be to compute the amount of 
J flowing into the BH and check for the conservation of 
J at the BH horizon (see e.g.. Sec. 4.2 of j90|). 



merger. The value of Jgum deviates from its initial value 
by about 2% after the merger at tret ^ 900M, but the 
deviation increases slowly and reaches ~ 7% at the end of 
simulation (tret = 2000M), indicating that the spurious 
loss of J continues after merger in case B4. 

The BH spin parameter for case B4 increases from 
aBu/MBu =0.75 to ~ 0.85 during the simulation (see 
Table [TTT|) . A secular decrease of total angular momen- 
tum as found in case B4 is commonly observed in rapidly 
spinning vacuum BH simulations (see e.g., [104j ), which 
find that spurious J loss is reduced with increasing reso- 



V. SUMMARY AND FUTURE WORK 

We present preliminary results from magnetized sim- 
ulations of BHNS binaries with BH:NS mass ratio 
3:1. We treat both initially nonspinning (cases A) and 
moderately-spinning (cases B) BHs. For those magnetic 
field configurations we consider, only initial NS magnetic 
fields with maximum (internal) strength of ~ lO^^G - 
corresponding to average magnetic to gas pressure ratio 
of 0.5% - significantly impact the inspiral and merger dy- 
namics. During merger, most of the magnetized NS mat- 
ter is captured by the BH. After disruption, the dynam- 



23 




tret/M f^^ [Hz] 



FIG. 22. Same as Fig.[2T]but for cases BO and B4. 



ics are followed for about 30-50 (Mo/1.4M0)ms. Only in 
the cases in which magnetic fields are strongest are mag- 
netic effects dynamically significant, increasing the final 
disk mass by up to a factor of two. The strong magnetic 
fields help push out the outer layers of the NS during tidal 
disruption, resulting in a gravitational wave mismatch of 
0.004 for the Advanced LIGO broadband configuration. 
It may be challenging for the upcoming generation of 
gravitational wave detectors to observe effects from such 
magnetic fields. Further studies with different field ge- 
ometries, black hole spins and higher resolution may be 
required to confirm this finding. 

Some SGRB models require a massive, hot, magne- 
tized disk around a BH with collimated magnetic fields to 
launch jets that generate 7-rays (see, e.g. [Hsl,!!^)- In our 
BHNS simulations, the remnant disk is hot (T ~ IMeV) 
and massive (Mdisk ~ O.O2M0 and ^ O.IMq for cases A 
and B, respectively), and possesses magnetic fields that 
are tightly wound. However, evidence for magnetic col- 
limation around the BH or magnetic field-induced tur- 
bulence in the disk is not observed. Future analyses will 
focus on mode growth studies and B-field decomposition 
in poloidal and toroidal components with respect to the 
centre of mass of the NS, similar to previous studies of 
single rotating or magnetized stars [l OSnlOTj . The lack 
of collimation may be due to the short disk evolution 
time before our simulation is terminated. The absence 
of turbulence in the disk may be due to insufficient res- 
olution in the disk, thereby suppressing instabilities like 
MRI. Therefore, future work will focus on evolving the 
disk at higher resolution, coupled with longer disk evolu- 
tions and different initial magnetic field configurations to 
more thoroughly assess the possibility of BHNS binaries 



FIG. 23. Gravitational- wave power spectrum for case AO com- 
puted as in Sec. IIIF of [5^. The solid curve shows the power 
spectrum of the numerical signal (dotted curve) stitched to 
the post-Newtonian waveform, including only the dominant 
(2, 2) and (2,-2) modes. The dashed curve shows the power 
spectrum expected from BHBH binaries with the same total 
mass and mass ratio as the BHNS, derived in [9^ from analy- 
sis of multi-orbit, non-precessing BHBH binaries. The heavy 
solid curve displays the effective strain of the Advanced LIGO 
detector, defined such that /iligo(/) = yfShXJ). The noise 
curve used here corresponds to the Advanced LIGO configu- 
ration ZERO_DET_HIGH_P [99]. To set physical units, a NS 
rest mass of Mo = 1.4M0 and a source distance of D=100Mpc 
is assumed. 



as short-hard GRB progenitors. 
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FIG. 25. Evolution of interior angular momentum Jint, angu- 
lar momentum carried off by GW A Jew, and BH's angular 
momentum Jbh for case A4. All quantities are normalized by 
the ADM angular momentum of the binary J. 



Appendix A: Hydrodynamic and MHD inequalities 

For a given set of "primitive" variables {p^^P^v^ ^B^) 
in the physical range (i.e. po > 0, P > and e > 0), 
the corresponding "conservative" variables {p^^f^Si^B^) 



must satisfy certain inequalities. In this appendix, we 
derive the inequalities and provide a practical recipe to 
impose these inequalities approximately to reduce inver- 
sion failures, which occur mainly in regions with very 
low density in the artificial "atmosphere" or inside the 
BH horizon where high accuracy is difficult to maintain 
but not crucial. Even when applying this recipe, inver- 
sion failures sometimes occur. In that case, we employ an 
alternative inversion scheme, described in Sec. IA4( that 
always works. Readers who are only interested in our 
recipe may jump directly to Sec. lA 31 and skip the rest of 
this appendix outlining the derivation. 



Since the inversion between B^ and B^ is trivial and 
does not involve other primitive variables, we will treat 
values of B^ as given and only consider the inversion 
from (p^^f^Si) to {pQ^P^v^). We assume that the EOS 
P = P{po^e) always gives P > whenever po > and 
e > 0. We also assume that the metric, lapse and shift are 
in the physical range. In particular, we require a > and 
jijk'^y > for any real vector k\ The requirement a > 
is always satisfied by our particular time slicing. How- 
ever, jij may lose positive-definiteness due to numberical 
error during the evolution, especially in the region deep 
inside the BH, near the "puncture". Before performing 
the inversion, we check if is positive definite by find- 
ing its eigenvalues. If jij is not positive definite, we reset 
jij i^^fij during the inversion, where fij is the 3D flat 
metric tensor. 
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1. Derivation of conservative variable inequalities: 
Pure hydrodynamic case 

In the absence of magnetic fields, the conservative vari- 
ables are given by 

P* = VlPolv (Al) 
Si = p^hui (A2) 
f = h{w- p*) + (/i - !)/>* - VtP, (A3) 



where 



= au 







(A4) 



and = 1 + e + P/po ^ 1- It follows from tx^tx^ = — 1 
that 



7^ ^Jl^Y^UiUj. 



(A5) 



Since Y^UiUj is positive definitive, 7^ > 1. We therefore 
conclude that 



p* > 



(A6) 



for po > and Wj G (—00,0x3). In addition, w = 7„/0* > 
, from which we obtain 



f>{h-l)p,-y/^P 
= V7bvpoe+{7v-l)P] >0. 



Hence we conclude that 



f > 0, 



(A7) 



(A8) 



which is Eq. (jT5|) . 

To derive the inequality (dH), Eqs. (|A2)) and (jA5]) are 
combined to yield 

52 = ^'^SiS, = {p.hfi^l - 1) = h\w^ - pi). (A9) 

A straightforward calculation yields 

(f + p,f =S^+ h'pl - y/^P{2hw - ^P) 

= ~S^ + {l.^f[pl{l + ef-P^] 

HVlPf ■ (AlO) 

Since there is no upper limit on 7^, the sum of the sec- 
ond and third terms is always positive if and only if the 
dominant energy condition - < Po(l + ~ holds. 
Hence we conclude that 



(All) 



if the dominant energy condition is satisfied. The in- 
equality (p!8|) is more stringent, and hence more useful. 
It can be derived using Eq. (jAlOp : 

f(f + 2p.) = 52 + (7„ V7)'(2pge + pie" - P') + P? . 

In this case, the sum of the second and third terms is 
always positive if and only if < 2pQe ^ p^e^ . Hence 
we have derived the inequality (Eq. [T8|): 



< f (f + 2p*) iff < 2p^e + pt,e\ 



(A12) 



Whether the condition P^ < 2pQe + pge^ is satisfied de- 
pends on the EOS. For the T-law EOS P = (E - l)poe, 
simple calulation gives 



P' -2pie-pie' =pie[T{T-2)e-2]. 



(A13) 



Hence the inequality S'^ < f(f + 2p*) holds iff r(r — 
2)e — 2 < 0. Restricting to the parameter space where 
the sound speed is subluminal, i.e. = VP/{pQh) < 1, 
we have r(r - 2)e < 1. Therefore, < f{f + 2p*) holds 
for the F-law EOS in regions where the sound speed is 
subluminal. For the F-law EOS it can be shown that 
< F — 1 holds for any nonnegative po and P. Thus, the 
sound speed will always be subluminal when 1 < F < 2. 
Thus, S'^ < f{f + 2p*) is satisfied for the F-law EOS, 
when 1< F < 2. 

We have just proved that the inequalities > 0, f > 
and 52 < f(f + 2p^) are necessary conditions for the 
primitive inversion to yield a physical solution for 1 < 
F < 2. We now want to prove that they are also the 
sufficient conditions for the F-law EOS with F > 1. For 
the F-law EOS, the enthalpy is related to the pressure P 

by 



TP 



(r-i)po' 



(AM) 



Combining the above equation with Eq. ()A3p yields 



Tw{f + p,)-{T-l)pl 

r«;2 _ (r _ i)p2 • 



(A15) 



It is useful to define a variable x = {w — p^)/f. It follows 
from Eqs. ()A15|) and (|A9]) that 



Ff(l — x){xf + p*) 
Tfx{fx + 2p*) + pI 



(AI6) 



and 



f{x) = fx{fx + 2p«) 



/l2 



0. 



(A17) 



These two equations can be combined to yield a quartic 
equation in x. However, it is easier to analyze the equa- 
tions in the present form. For any given p^ > 0, f > 
and S'^ < fir + 2p*), when x = 0, = 1 + Tf/p^ > 1 
and /(O) = -Syh'^ < 0; when x = 1, = 1 and 
/(I) = f{f-\-2p^) — S'^ >0. Hence the intermediate value 
theorem implies there exists xq G [0, 1] so that f{xo) = 0. 
The primitive variables are recovered from the following 
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expressions 



W TXq 

7^ = — = 1 + , 

Po 



h-1 



P = 



rf(i - x^)(x^f + p^) 
rfxo(fxo + 2p*) + p2 ' 

r-i 

— — Po(^- 1), 



(A18) 
(A19) 

(A20) 

(A21) 

(A22) 
(A23) 



Upon inspection, all the recovered primitive variables lie 
in the physically acceptable range for xq G [0, 1] and T > 
1. 

We therefore conclude that the inequalities p* > 0, 
f > and S'^ < f(f + 2p*) are both necessary and suf- 
ficient conditions for the inversion to yield physically- 
acceptable primitive variables for the F-law EOS with 
1 < r < 2. Since the T-law EOS with T = 2 is adopted 
in our simulations, we impose these inequalities in regions 
where there are no magnetic fields. 



2. Derivation of conservative variable inequalities: 
MHD case 



In the presence of magnetic fields, the conservative 
variables Si and f are given by 



(A24) 
(A25) 



Here Sf™'^ = p^hui and ffluid = h{w - p*) + {h — l)p* — 
-JjP, which are the same expressions as Eqs. (|A2|) and 
(|A3p . The variable = 7„ ^7 po remains unchanged and 
hence the inequahty > still holds in the presence of 
magnetic fields. 

It is convenient to introduce the following three quan- 
tities 



W = jyPoh = wh/y^ ■- 



ffluid 



{p.hf 



V7 



(A26) 

(A27) 
(A28) 

Following the algebra in Sec. 3.1 of J83, one can show 
that (c.f. Eqs. (26), (27) and (29) of [sf ) 

B'Si = (A29) 

{B^Sif{B^ + 2W) 



= {y/^fv\W + B^f 
V7 



r = rfl,a + Y(l + ^^)5^-^^, 



-,(A30) 
(A31) 



where B"^ = jijE^B^. It follows from Eq. (|X29|l and the 
Cauchy-Schwartz inequality that 

(B^Sif = {^'iBiS^"'^f 

< (7»JBi5j)(7»J5f"*<^5j'"*<^) 



= ^2^2 



fluid- 



Hence we obtain 



'S'fluid ^ {B^Sif, 



(A32) 



where = B^B. 

Using Eq. (|A30p . we can write 



and 



7(W^ + B2)2 



,2 {B^SiY{B' + 2W) 



(A33) 



Given values of Si and 5* , the only independent variable 
in the above equation is W . Straightforward calculation 
yields 



>0, 



(A35) 



where we have applied the Cauchy-Schwarz inequality 
{B'~Sif = (Y^BiSjf < {Y^BiBj){YiSiSj) = B^S\ 

Hence the maximum value of -Sl^-^ is achieved when 
ly ^ oo, which gives -Sl^-^ < S'^. The minimum value 
of -Sl^-^ is achieved when W = Wm^ where Wm is the 
minimum value of W for given values of p*. Si and 5*. 
Hence we obtain 



^ *^fluid ^ ^ •> 



where 



-2 ^ WlS^ + {WSif{B-^+2Wra) 

{W^ + B^f 



(A36) 



(A37) 



The upper and lower bounds of ffluid can be derived 
by first combining Eqs. ()A3ip and ()A33p : 



^^^'^-^ 2'' 2^7 (1^ + 52)2 • 



(A38) 



For fixed values of f. Si and 5% ffluid increases with 
increasing W . Therefore, we conclude that 



where 



~ ^ ~ ^ ~ a/7 o2 

Tm S Tfluid St —1^ , 



2 2^7 (W^™ + 52)2- 



(A39) 



(A40) 
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To calculate W^, we can combine the last equality in 
Eq. ()A26p and Eq. (|A34p to obtain an implicit relation 
of 1^ = W{h). Using V < 1, it is then straightforward 
to show that dW/dh > 0, and hence W = Wm when h is 
minimized. Given that /i > 1, we therefore have 



Wr.= 



■pi 



V7 



(A41) 



This equation can be combined with Eq. (jA37p . resulting 
in a quartic equation for Wm'- 

{lWl,-pl){Wra^B'f-WlS'^-{&~Sif{B'^2Wra) = 0, 

which may be solved analytically. Alternatively, 
Eqs. ()A37p and ()A4ip may be solved using an iterative 
scheme. We start by using Eq. ()A32p and choose an ini- 
tial guess Wm = + Next we compute 
the initial guess using Eq. (|A37j) . We then recompute 
Wm using Eq. (|A41[) and using Eq. (jA37j). We keep 
iterating until the values of and Wm converge. 

In the previous subsection, we proved that in the ab- 
sence of magnetic fields > 0, f > and S'^ < f(f+2p*) 
are the necessary and sufficient conditions for the inver- 
sion to produce the primitive variables in the physical 
range for the F-law EOS with 1 < F < 2. In the presence 
of magnetic fields, the necessary and sufficient conditions 
for the F-law EOS with 1 < F < 2 are replaced by the 
following inequalities: 

> 0, (A42) 
Tfluid > 0, (A43) 
Tfluid (Tfluid + 2p*) > S'l^id. (A44) 

Unfortunately, no simple, analytic expression for nec- 
essary and sufficient conditions between the conservatives 
Si^ f seems to exist in the presence of magnetic fields. 
However, necessary and sufficient conditions can be de- 
rived separately by combining Eqs. ()A42p - (|A44p . (|A36p 
and (|A39[) . The results are as follows. 

Necessary conditions for guaranteeing a physical solu- 
tion: 



p* > 0, 
- 2 



(A45) 
(A46) 



Sl<{f-^B^^ (f-^B'+2p^f. (A47) 

Sufficient conditions for guaranteeing a physical solu- 
tion: 

> 0, (A48) 
fm > 0, (A49) 
fm{rm^2p,)>S^. (A50) 

Both S'^ and fm are nonlinear functions of p*, 5^, and 
f. Unlike the pure hydrodynamic case, these inequal- 
ities are not trivial to impose strictly, so they are im- 
posed approximately as follows. First, a parameter 



is introduced, which determines whether the region under 
consideration is deep inside the BH horizon. For regions 
deep inside the BH horizon, defined by = > V^thr' 
the primary goal is to keep the evolution stable and pre- 
vent inaccurate data from leaking out of the BH horizon. 
We find that imposing the sufficient conditions ()A48p - 
(|A5Q[) approximately in this region is adequate (detailed 
recipe below). In regions where ?/^^ < V^^^r' So^l is 
to evolve the GRMHD equations as accurately as possi- 
ble, which means that imposing the sufficient conditions 
is not appropriate. We instead impose two of the nec- 
essary conditions ()A45p and ()A46p only. Since we do 
not strictly impose all the inequalities, inversion failures 
sometimes occur. Failures are fixed by replacing the f 
equation by the equation P = Pcoid(po)- We will demon- 
strate in Sec. IA4l that the set of equations — Jv\^ POi 



Eq. (|A24[) and P = Pcoid(po) always results in the prim- 
itive variables in the physical range. Our detailed recipe 
is described in the following subsection. 



3. Algorithm for Imposing MHD/Hydrodynamic 
Inequalities 

1. In any region, if < 0, set po = patm, P = ^atm, 
m = and recompute the conservative variables. 
If f < v/7 5V2, reset f = fatm + 3^2. 

2. In the region where 5* = 0, if f < 0, reset f = fatm- 
If S'^ > f{f + 2/)*), replace 



Si Si\ 



'f(f + 2p*) 



S2 



(A51) 



3. In the region where ip^ > ip^Yir (deep inside the BH 
horizon): First, estimate Wm and as follows 



W. 



mO 



1/2 



'^mO — 



{WmO + B^? 



1/2 



(A52) 
, (A53) 
(A54) 



Next, calculate Sm and fm from Eqs. ()A37P and 
(|A40|I . Then check if fm > fatm- If Tm < Tatm reset 
f (without changing Wm and Si) according to 



^B' 



Q2g2 _ (^B'SiY 
2^{Wm + B^?' 



(A55) 



Then check if S"^ < fm{fm + 2/9*). If not, reset Si 
(without changing fm and p*) according to 



Si ^ Si\ 



'fm{fm + 2/0*) 



52 



(A56) 
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These procedures do not strictly impose the suffi- 
cient conditions (|A49p and (|A5Qp , but can signifi- 
cantly reduce inversion failures. 

4. If the inversion still fails after going through all 
the above steps, replace the f equation ()A25p by 
P = Pcoid(Po) and perform the inversion. This 
procedure guarantees to produce the primitive vari- 
ables in the physical range, as will be shown in the 
next subsection. 



^coid(po) is continuous in [/?i,p2], which is true for the 
F-law EOS. This proves that the inversion always pro- 
duces po, as well as P in the physical range. The value 
of is then given by P^Ky/l Po)- For po ^ [pi,P2], 

we have 1 < 7^ < -^1 + S'^ / pi-, which is also in the 
physical range. Finally, the velocity is recovered from 
Eqs. (|A22]) and (|A23l) . 

Next consider the case ^ 0. Equations (jA26p and 
(IA341) yield 



4. Inversion using p*, Si and P = Pcoid(po) 

After imposing the conservative variables inequalities 
as described in the previous subsection, sometimes the 
conservatives^primitives variable inversion still fails. In 
this case, an alternative inversion is imposed, solving 
for po and Ui from the equations = s/jjvpo^ P = 
^coid(po) and Eq. ()A24p . We wih prove that this in- 
version always results in the primitive variables in the 
physical range. 

First consider the case 5* = for the F-law EOS with 
1 < F < 2. In this case, simple, analytic expressions for 
the necessary and sufficient conditions exist, are easy to 
implement (as described in Sees. [AT] fc lA 3[) . and guar- 
antee successful inversions. However, since the analysis 
is much simpler than the general case, it is instructive 
to study the alternative inversion scheme for this case 
ffist. It follows from = y/^Jvpo^ P = ^coid(po) and 
Eq. dAll) that 



Po 



p^ 
V7 



1-1/2 



P5^cold(P0) 



(A57) 



where /icoid = 1 + ecoid(po) + ^coid(po)/po is the specific 
enthalpy for the cold EOS. The above equation is an 
implicit equation for po- Define the function 



/(po) = Po 



p* 
V7 



1 



P*^cold(P0) 



-1/2 



(A58) 



and introduce two variables 

-1/2 

P* / . ^'M 

Pi 




P2 



p^ 



(A59) 



Since /icoid ^ 1, for given values of > and 5*^ G 
(—00,00), f{pi) < and f{p2) ^ 0. Hence there ex- 
ists po G [pi,P2] that satisfies Eq. (|A57[) provided that 



W 



q2 

O4 



^fluid + (P*^cold)^ 



fluid 



po 



(W + 52)2 



P* 



1-1/2 



fluid 



P5^cold(P0) 



(A60) 

(A61) 
(A62) 



where S^^^^ is regarded as an implicit function of Si, 
B\ p^ and po through Eqs. (jAGOp and ()A6ip . Hence 
Eq. ()A62p is an implicit equation for po- Next define 



/(po) = po 



p* 



02 

^fluid 



P5^cold(P0) 



-1/2 



(A63) 



where pi and p2 are as in Eq. ()A59p . Section [A3l proved 

that (B'Si)^ < Sl^.^ < S^ [Eqs. (lA32]) and (jA36l)]. It 
follows that /(pi) < and /(P2) > 0. Hence a solution 
exists for po G [pi,P2] as in the pure hydro case, pro- 
vided that both /icoid(po) and 5^|uid(Po) are continuous 
in [pi,p2]. The inversion thus produces po as well as P 
and in the physical range. Applying some algebraic 
manipulations to Eqs. ()A24p and (jSj, we recover the 4- 
velocity using the formula (c.f. Eq. (31) of (STf ) 



7vP*^cold(Po) 



Bi 



P*^cold(po) 



v/752 



and the 3- velocity from Eq. ()A23p . 

In conclusion, when using the equations p* = -s/y^yvpo^ 
P = Pcoid(po) and Eq. (|A24|) . the inversion will always 
produce the primitive variables in the physical range for 
any p* > and Si G (—00, 00). 
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